#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <sys/time.h>
#include <time.h>


#define EPSILON 1.0e-5f
#define REAL_T float
#define C 0.2f

/* Pour obtenir le temps "horloge" : */
double my_gettimeofday(){
	struct timeval tv;
	gettimeofday(&tv, NULL);
	return tv.tv_sec + (tv.tv_usec * 1.0e-6L);
}

static inline REAL_T calcul_iter_no_tests(REAL_T **grid[2], int grid_nb, int i, int j) {
	return ((1.0f - 4.0f * C) * grid[grid_nb][i][j]
	  + C * (grid[grid_nb][i][j-1])
	  + C * (grid[grid_nb][i][j+1])
	  + C * (grid[grid_nb][i-1][j]) 
	  + C * (grid[grid_nb][i+1][j]));
}

static inline REAL_T calcul_iter_tests (REAL_T **grid[2], int grid_nb, int i, int j, int h, int w) {
	return ((1.0f - 4.0f * C) * grid[grid_nb][i][j]
	  + C * (j>0   ? grid[grid_nb][i][j-1] : grid[grid_nb][i][j])
	  + C * (j<w-1 ? grid[grid_nb][i][j+1] : grid[grid_nb][i][j])
	  + C * (i>0   ? grid[grid_nb][i-1][j] : grid[grid_nb][i][j]) 
	  + C * (i<h-1 ? grid[grid_nb][i+1][j] : grid[grid_nb][i][j]));
}

int main(int argc, char **argv){
	int n = 0;
	int w = 0;
  	int h = 0;
  	int i = 0, j = 0;
  	int step = 0;

	int grid_nb = 0, inv_grid_nb=0;
  	REAL_T **grid[2] = {NULL, NULL};

	REAL_T max_error = 0.0;

	double t_debut = 0.0, t_fin = 0.0;
  	double t_total_debut = 0.0, t_total_fin = 0.0;
  

	/* Valeurs par defaut : */
  	w = 3;
  	h = 500;

  	/* Parametres : */
  	if( argc > 1) h    = atoi(argv[1]);
  	if( argc > 2) w    = atoi(argv[2]);

	/* Demarrage horloge : */
  	t_total_debut = my_gettimeofday();

  	/* Affichage informations diverses : */
  	printf("Dims : (%d, %d) \t Nb points : %i \t Taille memoire : %f Mo X2\n",
	 w, h, w*h, (float) w*h*sizeof(REAL_T) / (1024*1024));

  	/* Allocation memoire : */
  	for (grid_nb=0; grid_nb<=1; grid_nb++){
  	  	grid[grid_nb] = (REAL_T **) malloc(h*sizeof(REAL_T *));
  	  	if (grid[grid_nb] == NULL){
  	  		fprintf(stderr, "Echec allocation memoire pour 'grid[%d]'.\n", grid_nb);
  	  	}
		for (i=0; i<h; i++){
  	    	grid[grid_nb][i] = (REAL_T *) malloc(w*sizeof(REAL_T));
  	   	if (grid[grid_nb][i] == NULL){
				fprintf(stderr, "Echec allocation memoire pour 'grid[%d][%d]'.\n", grid_nb, i);
			}
      
			/* set to 0: */
      	memset(grid[grid_nb][i], 0, w*sizeof(REAL_T));
   	} /* for i */
  	} /* for grid_nb */

  	/* Initialisation par defaut de la grille 0 : */
  	grid_nb=0;
  	for (i=0; i<h; i++){
  		for (j=0; j<w; j++){
      	grid[grid_nb][i][j] = (REAL_T) 25.0;
    	}
  	}

  	/* Autres valeurs : 1 ligne de points a 100 toutes les 100 lignes */ 
  	step = 100;
  	for (i=step/2; i<h; i+=step){
  		for (j=0; j<w; j++){
    		if(j==0)
	    		printf("grid[grid_nb][%d][%d] = 100.0;\n", i, j);
      	grid[grid_nb][i][j] = 100.0;
    	} /* for j */
  	} /* for i */
  

  
   /* Test Exhaustif - Enregistrement des résultats */
	/*FILE * fi = fopen("results_seq_init", "w");
	 	 
	 for (i=0; i<h; i++){
	 	for(j=0;j<w;j++) {
	    	fprintf(fi, "%3.2f ", grid[grid_nb][i][j]);
   		}
   	
   		fprintf(fi, "\n");
    }
    
    fclose(fi);*/
  
  	/***** Calcul : *****/
  	do {
		inv_grid_nb = 1-grid_nb;

    	/* Debut horloge pour cette iteration : */
    	t_debut = my_gettimeofday();

    	/* mise a zero de grid[inv_grid_nb]: */
    	for (i=0; i<h; i++){
     		memset(grid[inv_grid_nb][i], 0, w*sizeof(REAL_T));    
    	} /* for i */

    	for (i=0; i<h; i++){
      
      	if(i==0 || i== h-1) {
				for (j=0; j<w; j++){
					grid[inv_grid_nb][i][j] = calcul_iter_tests(grid, grid_nb, i, j, h, w); 
	     	 	} /* for j */
   	   }
      	else
      	{
      		grid[inv_grid_nb][i][0] = calcul_iter_tests(grid, grid_nb, i, 0, h, w); 
      	
				for (j=1; j<w-1; j++){
					grid[inv_grid_nb][i][j] = calcul_iter_no_tests(grid, grid_nb, i, j); 
	     	 	} /* for j */
      	
      		grid[inv_grid_nb][i][w-1] = calcul_iter_tests(grid, grid_nb, i, w-1, h, w); 
      	}
		} /* for i */

		/* calcul de l'erreur */
   	max_error = 0.0f;
    	for (i=0; i<h; i++){
     		for (j=0; j<w; j++){
				REAL_T error = fabs(grid[grid_nb][i][j] - grid[inv_grid_nb][i][j]);
				if (error > max_error){ 
	  				max_error = error; 
				}
			} /* for j */
		} /* for i */

    	/* Arret horloge pour cette iteration : */
    	t_fin = my_gettimeofday();

    	/* Affichage pour cette iteration : */
    	printf("Iteration %d : \t delta = %.3e \t temps = %.1f s\n", 
		  	 ++n, max_error, t_fin - t_debut);
	
    	/* update 'grid_nb': */
    	grid_nb = inv_grid_nb;
	} while (max_error > EPSILON && n < 40);
  	/***** Fin du calcul *****/

	/* Test Exhaustif - Enregistrement des résultats */
	FILE * f = fopen("results_seq_opti", "w");
	 	 
	for (i=0; i<h; i++){
		for(j=0;j<w;j++) {
	   	fprintf(f, "%3.5f ", grid[inv_grid_nb][i][j]);
   	}
   	
   	fprintf(f, "\n");
	}
    
   fclose(f);

	/* Liberation memoire : */
  	for (grid_nb=0; grid_nb<=1; grid_nb++){
  		for (i=0; i<h; i++){
     		free(grid[grid_nb][i]);
    	}
  		free(grid[grid_nb]);
  	}

  	/* Arret horloge : */
  	t_total_fin = my_gettimeofday();  

  /* Affichage des resultats : */
  	printf("\nNombre d'iterations : %d\n", n); 
  	printf("Erreur (delta) = %.3e\n", max_error); 
	printf("Temps total : %.1f s\n", t_total_fin - t_total_debut); 

	return EXIT_SUCCESS;
}


